﻿Imports WINX32Lib 'Remember to add r3eference first (Add Reference option). For MS Visual Studio .NET 2003: before adding a reference to the project, corresponding assembly should be added to the Add Reference dialog box (http://msdn.microsoft.com/en-us/library/wkze6zky%28v=vs.71%29.aspx).
Imports System.IO
Imports System.Math

Module Analysis

    'this sub is defined to ovecome the VB.Net problem of coping arrays of different rank
    Declare Sub CopyMemory Lib "kernel32" Alias "RtlMoveMemory" _
        (ByVal Destination As Long, ByVal Source As Long, ByVal Length As Long)

    Sub Analyze(ByVal file_name As String, ByVal Frame_To_Analyze As Array, ByVal Total_Number_Of_Frames As Integer, _
             ByVal low As Short, ByVal high As Short, ByVal Min_Intensity As Short, ByVal low1 As Short, ByVal high1 As Short, _
             ByVal Detector_X_Dimention As Integer, ByVal Detector_Y_Dimention As Integer)
        Static Dim tic As String = 0
        Dim toc As String
        'tic = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60

        Static Dim Spectrum_File As New DocFile
        Static Dim Total_Number_Of_Good_Events_Per_Frame_File As New DocFile
        Static Dim Good_Events_File As New DocFile
        Static Dim Histogram_File As New DocFile
        Static Dim Multiple_Hits_File As New DocFile
        Static Dim Total_Number_Of_Multiple_Hits_File As New DocFile

        Static Dim Spectrum(,) As Single
        Static Dim Total_Number_Of_Good_Events(,) As Integer
        Static Dim Good_Events()(,) As Single 'Integer
        Static Dim Multiple_Hits()(,) As Integer
        Static Dim Total_Number_Of_Multiple_Hits(,) As Integer
        Static Dim Intensity_Histograme()() As Integer
        Static Dim Frame_Analized As Integer = 0

        Dim Max_Multiplicity As Integer = low \ (high - low)
        Dim Hist_Length As Integer = 100000 ' Min(Max_Multiplicity, 2) * high

        If Spectrum_File.GetParam(DM_CMD.DM_FILETITLE) Is Nothing Then
            Dim Spectrum_File_Path As String = Path.GetDirectoryName(file_name)
            Dim Spectrum_File_Name As String = Path.GetFileNameWithoutExtension(file_name)
            Total_Number_Of_Good_Events_Per_Frame_File.Open("", Total_Number_Of_Frames, 2, 1, dataType.X_LONG, Spectrum_File_Path & "\" & Spectrum_File_Name & "l.spe")
            Spectrum_File.Open("", Detector_X_Dimention, 2, 1, dataType.X_FLOAT, Spectrum_File_Path & "\" & Spectrum_File_Name & "s.spe")
            Good_Events_File.Open("", Detector_X_Dimention, Detector_Y_Dimention, 2, dataType.X_FLOAT, Spectrum_File_Path & "\" & Spectrum_File_Name & "o.spe")
            'Multiple_Hits_File.Open("", Detector_X_Dimention, Detector_Y_Dimention, Max_Multiplicity + 1, dataType.X_LONG, Spectrum_File_Path & "\" & Spectrum_File_Name & "mo.spe")
            Total_Number_Of_Multiple_Hits_File.Open("", Total_Number_Of_Frames, Max_Multiplicity + 1, 1, dataType.X_LONG, Spectrum_File_Path & "\" & Spectrum_File_Name & "ml.spe")
            Histogram_File.Open("", Hist_Length/10, 1, 3, dataType.X_LONG, Spectrum_File_Path & "\" & Spectrum_File_Name & "h.spe")


            ReDim Spectrum(Detector_X_Dimention - 1, 1)
            ReDim Multiple_Hits(Max_Multiplicity)
            For i As Integer = 0 To Max_Multiplicity
                ReDim Multiple_Hits(i)(Detector_X_Dimention - 1, Detector_Y_Dimention - 1)
            Next
            ReDim Preserve Total_Number_Of_Multiple_Hits(Total_Number_Of_Frames - 1, Max_Multiplicity)
            ReDim Preserve Good_Events(1)
            ReDim Preserve Good_Events(0)(Detector_X_Dimention - 1, Detector_Y_Dimention - 1)
            ReDim Preserve Good_Events(1)(Detector_X_Dimention - 1, Detector_Y_Dimention - 1)
            ReDim Preserve Total_Number_Of_Good_Events(Total_Number_Of_Frames - 1, 1)
            ReDim Preserve Intensity_Histograme(2)
            ReDim Preserve Intensity_Histograme(0)(Hist_Length) '(Max_Multiplicity * high)
            ReDim Preserve Intensity_Histograme(1)(Hist_Length) '(Max_Multiplicity * high)
            ReDim Preserve Intensity_Histograme(2)(Hist_Length) '(Max_Multiplicity * high)
        End If


        If Frame_To_Analyze Is Nothing Then
            Frame_Analized = 0
            Spectrum = Nothing
            Good_Events = Nothing
            Multiple_Hits = Nothing
            Intensity_Histograme = Nothing
            Exit Sub
        Else
            Frame_Analized = Frame_Analized + 1


            Dim toc2, toc3 As String
            toc3 = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60 - tic
            tic = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60
            '%%%%%%%%%%%%%%%%%%%%%%
            'Strange behaviour - the function call takes more than 1s
            Analyze_Frame_Grains(Frame_To_Analyze, Detector_X_Dimention - 1, Detector_Y_Dimention - 1, Frame_Analized - 1, _
                          low, high, Min_Intensity, low1, high1, _
                          Spectrum, Good_Events, _
                          Total_Number_Of_Good_Events, Intensity_Histograme, _
                          Multiple_Hits, Total_Number_Of_Multiple_Hits, toc2)
            toc = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60 - tic
            tic = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60
            '%%%%%%%%%%%%%%%%%%


            Spectrum_File.PutFrame(1, Spectrum)
            Good_Events_File.PutFrame(1, Good_Events(0))
            Good_Events_File.PutFrame(2, Good_Events(1))
            Total_Number_Of_Good_Events_Per_Frame_File.PutFrame(1, Total_Number_Of_Good_Events)
            Total_Number_Of_Multiple_Hits_File.PutFrame(1, Total_Number_Of_Multiple_Hits)

            'For i As Integer = 0 To Max_Multiplicity
            '    'Dim tem_Multiple_Hits(Multiple_Hits.GetLength(0) - 1, Multiple_Hits.GetLength(1) - 1) As Array
            '    'Array.Copy(Multiple_Hits, i * tem_Multiple_Hits.Length, tem_Multiple_Hits, 0, tem_Multiple_Hits.Length)
            '    Multiple_Hits_File.PutFrame(i + 1, Multiple_Hits(i))
            'Next

            Dim Hist_To_Plot(,) As Integer
            If Not Frame_Analized = Total_Number_Of_Frames Then
                ReDim Hist_To_Plot(Hist_Length / 10 - 1, 0)
                ''we dont need to show the 0 intensity in the histograme - we start at 1
                'Buffer.BlockCopy(Intensity_Histograme(0), 1 * 4, Hist_To_Plot, 0, Hist_Length * 4)
                'Buffer.BlockCopy(Intensity_Histograme(1), 1 * 4, Hist_To_Plot, Hist_Length * 4, Hist_Length * 4)
                'Buffer.BlockCopy(Intensity_Histograme(2), 1 * 4, Hist_To_Plot, Hist_Length * 4 * 2, Hist_Length * 4)
                'Histogram_File.PutFrame(1, Hist_To_Plot)
                Buffer.BlockCopy(Intensity_Histograme(0), 4, Hist_To_Plot, 0, Hist_Length / 10 * 4)
                Histogram_File.PutFrame(1, Hist_To_Plot)
                Buffer.BlockCopy(Intensity_Histograme(1), 4, Hist_To_Plot, 0, Hist_Length / 10 * 4)
                Histogram_File.PutFrame(2, Hist_To_Plot)
                Buffer.BlockCopy(Intensity_Histograme(2), 4, Hist_To_Plot, 0, Hist_Length / 10 * 4)
                Histogram_File.PutFrame(3, Hist_To_Plot)
            Else
                Dim Final_Length As Integer = 10000 ' Array.FindLastIndex(Intensity_Histograme(0), AddressOf isNotZero)
                ReDim Hist_To_Plot(Final_Length - 1, 0)
                Dim Spectrum_File_Path As String = Path.GetDirectoryName(file_name)
                Dim Spectrum_File_Name As String = Path.GetFileNameWithoutExtension(file_name)
                Histogram_File.Open("", Final_Length, 1, 3, dataType.X_LONG, Spectrum_File_Path & "\" & Spectrum_File_Name & "h.spe")
                Buffer.BlockCopy(Intensity_Histograme(0), 4, Hist_To_Plot, 0, Final_Length * 4)
                Histogram_File.PutFrame(1, Hist_To_Plot)
                Buffer.BlockCopy(Intensity_Histograme(1), 4, Hist_To_Plot, 0, Final_Length * 4)
                Histogram_File.PutFrame(2, Hist_To_Plot)
                Buffer.BlockCopy(Intensity_Histograme(2), 4, Hist_To_Plot, 0, Final_Length * 4)
                Histogram_File.PutFrame(3, Hist_To_Plot)
            End If
            toc = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60 - tic
        End If


    End Sub

    Function isNotZero(ByVal a As Integer) As Boolean
        Return Not a.Equals(0)
    End Function

    Sub Clear_DocFile(ByRef Doc_File As DocFile)
        For i As Integer = 1 To Doc_File.GetParam(DM_CMD.DM_NUMFRAMES)
            Dim Temporary_Frame As Array
            Doc_File.GetFrame(i, Temporary_Frame)
            Array.Clear(Temporary_Frame, 0, Temporary_Frame.Length)
            Doc_File.PutFrame(i, Temporary_Frame)
        Next
    End Sub

    Sub Analyze_Frame(ByRef Frame_To_Analyze As Array, ByRef Frame_DimX As Integer, ByRef Frame_DimY As Integer, _
                     ByRef Frame_Index As Integer, _
                 ByVal low As Short, ByVal high As Short, ByVal Min_Intensity As Short, ByVal low1 As Short, ByVal high1 As Short, _
                 Optional ByRef Spectrum(,) As Integer = Nothing, _
                 Optional ByRef Good_Events_Sum()(,) As Integer = Nothing, _
                 Optional ByRef Total_Number_Of_Good_Events(,) As Integer = Nothing, _
                 Optional ByRef Intensity_Histograme()() As Integer = Nothing, _
                 Optional ByRef Multiple_Hits()(,) As Integer = Nothing, _
                 Optional ByRef Total_Number_Of_Multiple_Hits(,) As Integer = Nothing, Optional ByRef toc2 As String = "")

        Dim tic As String
        tic = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60

        'to speed up we memorize the following:
   


        If Not Frame_To_Analyze.Rank = 2 Then
            MsgBox("The frame to analize hase a dimention of " & Frame_To_Analyze.Rank)
        End If

        Dim Max_Multiplicity As Integer = low \ (high - low)
        Static Dim Max_Intensity As Integer = Intensity_Histograme(0).GetLength(0) - 1

        Dim Good_Events(,,) As Integer
        ReDim Good_Events(Frame_DimX, Frame_DimY, 1)

        Dim number_of_points_for_neighbour_analize As New Integer
        Dim neighbour_analize_index(,) As Integer
        ReDim neighbour_analize_index((Frame_DimY + 1) * (Frame_DimX + 1) - 1, 1)

        'to speed up the code we introduce Return_To_Neighbours Array for pixels with position i=0 or j=0
        Dim Return_To_Neighbours() As Boolean
        ReDim Return_To_Neighbours(Frame_DimY)

        'After having alook into jakubs code I decided to change the sope of the loop not to tuch the borders
        For i As Integer = 1 To Frame_DimX - 1
            For j As Integer = 1 To Frame_DimY - 1
                'we skip most of the loop if the intensity is to low (most frequent case) 
                If Frame_To_Analyze(i, j) < Min_Intensity Then
to_execute_if_the_intensity_is_too_low:
                    If Frame_To_Analyze(i, j) > 0 Then
                        Intensity_Histograme(0)(Frame_To_Analyze(i, j)) += 1
                        Intensity_Histograme(1)(Frame_To_Analyze(i, j)) += 1
                        'Else
                        '    Frame_To_Analyze(i, j) = 0
                    End If
                Else
to_execute_if_the_intensity_is_abowe_treshhold:
                    'the Intensity_Histograme array can grow as we don't know the maximum intensity value
                    If Max_Intensity < Frame_To_Analyze(i, j) Then
                        Max_Intensity = Frame_To_Analyze(i, j) * 2
                        ReDim Preserve Intensity_Histograme(0)(Max_Intensity)
                        ReDim Preserve Intensity_Histograme(1)(Max_Intensity)
                    End If
                    Intensity_Histograme(0)(Frame_To_Analyze(i, j)) += 1
                    Intensity_Histograme(1)(Frame_To_Analyze(i, j)) += 1
Start_of_the_multiple_hits_analysis:
                    'here we check the multiplicity of the hit
                    Dim Lower_Lewel_Multiplicity As Integer = Frame_To_Analyze(i, j) \ low
                    If Lower_Lewel_Multiplicity > 0 Then 'one should be creful to satisfy the condition : high < 2*low
                        Dim Upper_Lewel_Multiplicity As Integer = Ceiling(Frame_To_Analyze(i, j) / high)
                        If Lower_Lewel_Multiplicity >= Upper_Lewel_Multiplicity Then
                            'Good_Events(i, j, 0) = Lower_Lewel_Multiplicity
                            Spectrum(i, 0) += Lower_Lewel_Multiplicity
                            Good_Events_Sum(0)(i, j) += Lower_Lewel_Multiplicity
                            Total_Number_Of_Good_Events(Frame_Index, 0) += Lower_Lewel_Multiplicity
                            If Lower_Lewel_Multiplicity = Upper_Lewel_Multiplicity Then
                                Multiple_Hits(Lower_Lewel_Multiplicity)(i, j) += 1
                                Total_Number_Of_Multiple_Hits(Frame_Index, Lower_Lewel_Multiplicity) += 1
                            Else
                                'Good_Events(i, j, 1) = Upper_Lewel_Multiplicity - Lower_Lewel_Multiplicity
                                Spectrum(i, 1) += Upper_Lewel_Multiplicity - Lower_Lewel_Multiplicity
                                Good_Events_Sum(1)(i, j) += Upper_Lewel_Multiplicity - Lower_Lewel_Multiplicity
                                Total_Number_Of_Good_Events(Frame_Index, 1) += Good_Events(i, j, 1)
                            End If
                        End If
                        'for aprioprate neighbours analysis we need to clear out the  pixels with the intensity higher than low
                        Frame_To_Analyze(i, j) = 0
                    Else
                        neighbour_analize_index(number_of_points_for_neighbour_analize, 0) = i
                        neighbour_analize_index(number_of_points_for_neighbour_analize, 1) = j
                        number_of_points_for_neighbour_analize += 1
                    End If
                    'Spectrum(i, 0) = Spectrum(i, 0) + Good_Events(i, j, 0)
                    'Spectrum(i, 1) = Spectrum(i, 1) + Good_Events(i, j, 1)
                    'Good_Events_Sum(0)(i, j) = Good_Events_Sum(0)(i, j) + Good_Events(i, j, 0)
                    'Good_Events_Sum(1)(i, j) = Good_Events_Sum(1)(i, j) + Good_Events(i, j, 1)
                End If
            Next
        Next

        'neghbours_analysis
        Dim Neighbours_Sum As New Integer
        For k As Integer = 0 To number_of_points_for_neighbour_analize - 1
            '    If Not 0 = Frame_To_Analyze(neighbour_analize_index(k, 0), neighbour_analize_index(k, 1)) Then
            Neighbours_Sum = Frame_To_Analyze(neighbour_analize_index(k, 0) - 1, neighbour_analize_index(k, 1) - 1) + _
                             Frame_To_Analyze(neighbour_analize_index(k, 0) - 1, neighbour_analize_index(k, 1)) + _
                             Frame_To_Analyze(neighbour_analize_index(k, 0) - 1, neighbour_analize_index(k, 1) + 1) + _
                             Frame_To_Analyze(neighbour_analize_index(k, 0), neighbour_analize_index(k, 1) - 1) + _
                             Frame_To_Analyze(neighbour_analize_index(k, 0), neighbour_analize_index(k, 1)) + _
                             Frame_To_Analyze(neighbour_analize_index(k, 0), neighbour_analize_index(k, 1) + 1) + _
                             Frame_To_Analyze(neighbour_analize_index(k, 0) + 1, neighbour_analize_index(k, 1) - 1) + _
                             Frame_To_Analyze(neighbour_analize_index(k, 0) + 1, neighbour_analize_index(k, 1)) + _
                             Frame_To_Analyze(neighbour_analize_index(k, 0) + 1, neighbour_analize_index(k, 1) + 1)
            If (Neighbours_Sum < high1) AndAlso (Neighbours_Sum > low1) Then
                Intensity_Histograme(1)(Neighbours_Sum) += 1
                Intensity_Histograme(2)(Neighbours_Sum) += 1
                For i = neighbour_analize_index(k, 0) - 1 To neighbour_analize_index(k, 0) + 1
                    For j = neighbour_analize_index(k, 1) - 1 To neighbour_analize_index(k, 1) + 1
                        'to speed up we omit the pixels with the intensity<=0
                        If Not Frame_To_Analyze(i, j) <= 0 Then
                            Intensity_Histograme(1)(Frame_To_Analyze(i, j)) -= 1
                            Frame_To_Analyze(i, j) = 0
                        End If
                    Next
                Next
                Multiple_Hits(0)(neighbour_analize_index(k, 0), neighbour_analize_index(k, 1)) += 1
                Total_Number_Of_Multiple_Hits(Frame_Index, 0) += 1
                'Good_Events(i, j, 0) = Good_Events(i, j, 0) += 1 ' not needed
                Total_Number_Of_Good_Events(Frame_Index, 0) += 1
                Spectrum(neighbour_analize_index(k, 0), 0) += 1
                Good_Events_Sum(0)(neighbour_analize_index(k, 0), neighbour_analize_index(k, 1)) += 1
            End If
        Next
        toc2 = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60 - tic
    End Sub



End Module
